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We introduce the study of dynamical quantum noise in Bose-Einstein condensates through nu- 
merical simulation of stochastic partial differential equations obtained using phase space represen- 
tations. We derive evolution equations for a single trapped condensate in both the positive-P and 

■ Wi gner representations, and perform simulations to compare the predictions of the two methods. 
The positive-P approach is found to be highly susceptible to the stability problems that have been 
observed in other strongly nonlinear, weakly damped systems. Using the Wigner representation, we 

— ^ examine the evolution of several quantities of interest using from a variety of choices of initial state 

Q . for the condensate, and compare results to those for single-mode models. 

<^ • I. INTRODUCTION 

S: 

' A key focus of the explosion of interest in the new dilute atomic gas Bose-Einstein condensates [Q-H has been the 
(— I , study of the time-evolution of condensates from some initial state. Amongst many works, there have been theoretical 
Q ' investigations of the way condensates react to a range of perturbations, such as "shaking" the trap to excite sound 
O ] waves removing a potential barrier to allow two condensates to interfere applying electromagnetic fields 

to transfer condensate population into other possibly untrapped states [p[-p^, or "stirring" a condensate to excite 

■ vortices jlj] . All but the last of these effects have already been demonstrated experimentally. A common element in 
^ . the theoretical works on these topics is the description of the condensate using the time-dependent Gross-Pitaevskii 

0^ ■ equation (GPE) or coupled GPEs (or their approximate hydrodynamic versions). The GPE can be derived as an 
\l , equation for the condensate amplitude assuming the condensate state is a multi-mode coherent state (on the concept 

■ of coherent states see |l^ , p^ ). Hence an implicit assumption underlying these approaches is that the condensate 
is adequately described as a coherent state. However, a number of experiments are now exploring issues such as 
coherence [|6[0 and the diffusion of relative phase between two condensates These concepts are familiar 

_ from optical systems, but additional factors arise in condensate physics such as the dispersion associated with the 
^ . nonzero atomic mass, and especially the effects of atomic interactions. In particular, several early models for atom 
d ' lasers |l9|-^ suggest that coherence properties may be strongly influenced by the nonlinear interactions. Moreover, 
H one of the principal themes of quantum optics is the idea that processing of quantum noise by nonlinearities leads to 



I interesting statistical properties |14 15[|. It thus becomes important to consider the nature of the condensate beyond 
' ^ , the coherent state, and in particular the influence of quantum noise on the coherence of the condensate and any 
eventual "atom laser" . For example, we may ask how different loss profiles for an output coupler might affect the 
noise statistics of an atom laser? 

In fact, there have been a number of studies into intrinsically quantum dynamical effects, in particular the collapse 
and revival of the relative phase between two coupled condensates and the robustness of such effects against 

environmental decoherence |^,^. However, these studies have been restricted to one or two modes and assume the 
' condensate wave function is independent of the number of atoms. For large condensates, these approximations are 
5^ , not necessarily valid and a many-mode approach is required. A single mode model may give an estimate for the phase 
diffusion time for example, but can never describe spatial coherence properties or the role of local density and 
phase fluctuations. Hence, there is a need for techniques to treat the quantum dynamics of the condensate using a 
fully spatially dependent field rather than a few mode approach. 
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Closely related issues are well-known in the field of quantum optical solitons and nonlinear quantum optics in 
general. In that situation, the propagation of the optical field is governed by a quantum nonlinear Schrodinger 
equation that includes the effects of fiber dispersion and the Kerr nonlinearity of the medium ||2^,^ . Such a model 
leads to the prediction that a soliton pulse injected into the fiber experiences squeezing (for general discussions of 
,) in both the electric field amplitude |^,^ and the photon number |^,|l|. Both types of 



squeezmg see 

squeezing have now been observed [[}2|,|33[ . While estimates for the squeezing have also been obtained from single 
mode models of the Kerr nonlinearity p5| , the presence of fiber dispersion means that accurate results can only be 
obtained from a multi-mode treatment of the full quantum field. The Heisenberg field operator describing a BEC 
confined in a one-dimensional potential well obeys a strikingly similar equation to that of the fiber soliton system, 
differing only by the addition of the trapping potential and the interpretation of the dispersive term which in a 
condensate represents the kinetic energy. Thus given that multi-mode models have proved essential for the accurate 
prediction of quantum soliton properties, it is reasonable to assume the same may hold true in Bose condensates. In 
fact, we see below that the nonlinearity occurring in the condensate problem is typically far larger than for the soliton 
case, and thus the role of quantum noise should be more important. 

While we thus have strong motives for seeking the complete evolution of the field operator, such a calculation is 
at first sight a formidable task, if for no other reason than that the Hilbert space for the system is truly vast. The 
numerical calculation of the evolution of just a single operator with significant excitation requires a large basis. The 
problem of the field is far worse. Nevertheless, in this paper we intend to demonstrate how techniques of quantum 
optics may be used to provide a complete description of the condensate field operator such that we can calculate 
virtually any desired quantum expectation value. The key to our approach is the representation of the density 
operator using phase-space quasi-probability functions. These functions then lead naturally to a description in terms 
of classical fields which are subject to evolution equations similar to the semi-classical time-dependent Gross Pitaevskii 
equation (GPE) satisfied by the mean field but with the crucial addition of stochastic driving terms. These terms 
do not correspond to any physical noise sources, but are defined in such a way as to recapture exactly the operator 
character of the fully quantum-mechanical field. In particular we use two representations: the positive-P function and 
the well-known Wigner distribution. While we can perform exact calculations in the positive-P representation, we find 
the system rapidly succumbs to the instabilities that have been observed previously for that representation p^ , ^ . 
Therefore we also consider an approximate but robust method using the Wigner representation. We are then able to 
extract a large range of interesting averages. 

It is worth noting that in terms more familiar in conventional quantum field theory, the stochastic techniques 
we present in this paper constitute a method of numerically evaluating path-integral representations of quantum 
field averages (see Plimak et al js^-Q). The phase-space of the classical fields is in these terms the space of the 
Feynman paths, while the phase-space quasi-probability functions are measures over the respective path integrals. 
These measures are constructively characterised by the corresponding stochastic evolution equations, which allows 
the path integrals to be calculated. This point highlights the significance of the positive-P representation. Providing 
certain boundary conditions are satisfied ||lj,|l^ , the positive-P is an exact method for propagating the field in real (as 
opposed to imaginary,) time. While it may blow up in certain cases as time progresses, all other exact methods fail at 
all times! Direct integration of the Feynman path integrals for example, is numerically useless due to the oscillatory 
phase factors. 

The stochastic techniques are well-known in quantum optics, and may be familiar to readers with a background in 
that subject. However, many of these ideas may be new to those who have come to Bose-Einstein condensation from 
other disciplines. Therefore we have taken a pedagogic approach in deriving the fundamental stochastic equations. 
Before treating the full quantum field problem, we review the techniques in the single-mode approximation for which 
the Hamiltonian corresponds to the anharmonic oscillator. The analysis for the complete field then follows in a natural 
way. Readers already well-versed in the stochastic approach to quantum dynamics may wish to pass over this earlier 



background in section [IC. 

The paper is structured as follows. In section |l| we provide a detailed demonstration of the techniques f or p rop- 



II B, we 



agation of quantum fields using phase space representations. Af ter s tating the complete problem in section 
simplify to the corresponding one- mode Hamiltonian in section [I C and derive the equivalent sto chas tic equations 
in the positive-P representation in detail. We generalize this approach to the full field in section IIP , and find an 
approximate but more stable method using the Wigner function in section HE. There is considerable freedom in the 



choice of initial states for our simulations. We discuss these issues and present some natural choices in section IH 



Our numerical results illustrating some of the possibilities of the stochastic approach are given in a range of section IV 
before we conclude. 
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II. TECHNIQUES FOR PROPAGATION OF QUANTUM FIELDS 



A. General Ideas 



There are a number of well-established techniques in quantum optics for the propagation of a complete quantum 
field. Typically, these ideas involve a generalization of standard procedures for finding the time evolution of averages 
in a system with a single mode or a few modes p^ , |l5| . In summary, the procedure is as follows. The system 
density operator is expressed in the coherent state basis using a quasi-probability function such as the P, Wigner 
or positive-P distributions. The master equation describing the evolution of the density operator is converted to an 



equivalent partial differential equation (PDE) for the distribution. If certain conditions are satisfied 15|, the PDE 
may then be converted to a set of classical stochastic ordinary differential equations (Langevin equationsjwhich yield 
quantum expectation values as ensemble averages of moments of the phase space variables. This procedure has the 
significant advantage of providing a natural numerical implementation in which we calculate the evolution of a small 
number of phase space variables rather than that of a very large number of variables describing the density matrix 
or a distribution function on a large complex grid. This advantage becomes essential for systems of several modes for 
which the Hilbert space is so large that a direct numerical simulation would be impossible. Now when we consider a 
quantum field, unless we are fortunate enough to have an analytic solution, the problem must be treated numerically 
as a system with a large but finite number of modes and the associated Hilbert space is truly vast. A thousand atoms 
with just one hundred modes for example, occupy a Hilbert space of dimension 100^000. The stochastic treatment 
is now vital. One uses essentially the same procedure but works with functional PDE's, and thus obtains stochastic 
equations for a classical field ||2^,^,^,Q . In the following sections we provide a detailed derivation for the case of a 
Bose condensate in one dimension. 



B. Trapped Bose-Einstein condensates 

We model a one-dimensional system by assuming a highly anisotropic harmonic trap with the longitudinal and 
radial trap frequencies {lOz and LOr respectively ) satisfying A = ujz/uJr 1- Below, we use parameters corresponding 
to a cigar-shaped trap such as that in Rcfs. pUl^]. With strong radial confinement, we assume the nonlinearity plays 
a negligible role in the radial direction. The field operator is then assumed to factorize with its transverse dependence 
completely described by a coherent state occupation of the lowest mode of the trap. So the Heisenberg picture boson 
field operator has the form 

n^) = {^) exp -^U(.,t). (1) 



Adopting harmonic oscillator units in the axial direction with oq = y^h/{mujz), t — ujzt, x — z/aa and ^l;{x,T) = 
y/ao4>{z,t), the one-dimensional second-quantized Hamiltonian is 

/OC /•OO 
dxtP^ICi^ + - dxi^^^i^, (2) 
CO ^ J — oo 

where JC is the hnear operator 

^ = -^d^ + k-^' 

fi is the scaled chemical potential, and T = 2a/(Aao) is the scaled nonlinear constant with a the s-wave scattering 
length. Our ultimate aim is to calculate (multi-time) averages of ip under the evolution induced by the Hamiltonian (H). 
More generally, the system may include damping in the form of couplings to atom reservoirs. In this case the system 
is described by a density matrix satisfying a master equation 

where the Liouvillian £ is a superoperator that acts to the right in the fashion 

Cp = -z[H, p] + Y, y(20,p6] - d]d,p - pd]6,), (5) 

j 
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and the Oj are operators describing the bath coupUngs with strengths Kj. While we do not include damping in the 
present work, it is convenient to work in a density matrix formalism. We now develop stochastic descriptions of the 
dynamics in both the P and Wigner representations, explaining the method for the P representation in detail. 



C. One mode problem 

1. P representation 

To illustrate the ideas underlying a phase-space approach without the notational baggage of the full many-mode 
problem, we begin by treating the single mode limit of Eq. (H). The atomic field is assumed to be described by a 
single mode operator a{t) with an associated mode function -^gp determined by the solution to the time-independent 
GPE. The single mode Hamiltonian is 

IP^^'^ =ua)a+^a)a)aa, (6) 

with uj — da; i/'Sp^V'gp ^^nd x = F da; |?/'gp|'*- The first step in our attempt to obtain a stochastic description 
is to express the density matrix in a diagonal coherent state basis using the Glauber- Sudarshan P function p^ , p^ : 

p = j (fa\a){a\P[a), (7) 

where |a) is a coherent state with the c-number complex amplitude a. It is tempting to consider P{a) as a probability 
distribution for the density matrix over the coherent states. However, while P{a) is real, for non-classical states it 
may be highly singular and/or take on negative values ||l^,|l5|. It is proved in Ref. |0 that there is a unique P 
function for every density matrix. The quantum averages of interest are found as moments of the P distribution 
which correspond to norman?/-ordered expectation values: 

(at^a") = Tr{at™a"p} = ( d? a a*"" a" P {a) , (8) 



for integer m,n > 0. Arbitrarily-ordered averages can always be found by first rewriting them in terms of normally- 
ordered quantities. 

We now need an equation for the time-evolution of the P 
function. Using |a) = exp(— |q;P/2) exp(aa^) exp(— a+a) |0) and the definition of the P function it is not hard to 
demonstrate the operator correspondences: 

txp ^ aP{a), a) p ^ — 

pa ^ (a — ^ I Pia), pa)^a^P{a). (9) 

\ oa+ J 

We have introduced the unusual notation a+, which for the moment is to be read as the ordinary complex conjugate 
a*. Substituting these correspondences in the master equation p = —i[H^^\p] we obtain a Fokker-Planck equation 
for the time-evolution of P{a) 

dP ( d d 19^ 19^ 1 

Note that despite the appearance of this equation, a and are not to be treated as independent variables as they 
are complex conjugates jl4|Jl5[| . 

Equation ( p^ ) is exact and completely equivalent to the master equation. Our motivation in obtaining it is based 
on the fact that any Fokker-Planck equation with a positive-definite diffusion matrix may be exactly rewritten in the 
language of stochastic differential or Langevin equations To be precise, consider a Fokker-Planck equation of the 
form 

f - E iM-. 0^ 4 £ il-ilr I'*'- ')>^"<- 
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in which the diffusion matrix D = BB^ is clearly positive definite. Then a third equivalent description is given by 
the system of stochastic equations 

^=A(x,t)+B(x,t)E(t), (12) 



where the real noise sources Ej{t) have zero mean and satisfy Ej{t)Ek{t') ~ SjkS{t — t'). These equations (and all 
other stochastic equations in this paper) are to be interpreted in the Ito approach to stochastic calculus. A complete 
discussion of the techniques of stochastic calculus and the connection between the Fokker- Planck and Langevin 
descriptions is provided by Gardiner |4l| . By making such a transformation we would apparently have achieved our 
aim of a stochastic description of the quantum dynamics, and could calculate expectation values by taking ensemble 
averages of moments of the phase space variables xj. 



2. Positive-P representation 

Unfortunately, the diffusion matrix in Eq. (|l^) is clearly not positive-definite and the preceding equivalence does 
not apply. However, Drummond and Gardiner have shown that in such cases, the situation may be rescued 
by introducing the "Positive-P" function which represents the density matrix as an integral over two independent 
variables 

p = |dW/3M^p(«,;3). (13) 

It can be shown that with this definition the positive-P function is analytic and can be chosen positive for any 
density matrix The crucial step comes here. Referring to Eq. ([l0|), we now consider a and as independent 

quantities and by making the identification P — , we may read Eq. (|l^) as the Fokker-Planck equation in the 
positive-P representation corresponding to the original master equation. Moreover, by writing P(a, /?) as a function 
of four real variables rather than two complex variables, one may show that the resulting 4x4 diffusion matrix is 
always positive-definite Thus in the positive-P representation, it is always possible to derive an equivalent 

Langevin equation description using Eq. (p^. In the present case, Eq. ( p^ leads to the stochastic system 

*^ ^ + Xa^a^) + \/ixarii{t)^ 
^cfa ^ _(^^^+ y^aa'^'^) + ^/-ixa'^mit)^ (14) 

CET 

where rii{t) for i ~ 1,2 are delta-correlated in time with zero- mean. Note that a and a"*" experience different noise 
sources, so that even if they are conjugate at r = they do not remain so. 

We emphasize that Eqs. (^^ are completely equivalent to the master equation. Any expectation value that can 
be found from the density operator may equally be found by ensemble averaging over many trajectories using the 
stochastic equations. Being c-number variables, a and do not satisfy the commutation properties of the operators a 
and a^. Nevertheless, through the inclusion of the noise sources r/i and the insistence on normal ordering when taking 
averages, they still account for the complete quantum dynamics. This equivalence of the stochastic and operator 
approaches has been demonstrated explicitly in the context of optical fiber solitons in a recent paper of Fini et al p3| . 
It is important to note that the r]i do not correspond to any physical noise sources, but are included only to recapture 
the commutation relations of the operators. In this sense they are quite distinct from the operator valued noise sources 
that appear in "quantum Langevin equations" |l^jl5[| . In fact in the positive-P representation, Plimak et al [^sj have 
pointed out that there is some freedom in the precise form of the noise terms, which can be exploited to improve 
convergence properties dramatically. We also mention a well-known difficulty with the positive-P representation. 
The independence of the noise sources driving a and can in some cases lead to wild trajectories that prevent 
convergence of the ensemble averages |34|j35| , ^^ . This is indeed true in the present case and we consider this problem 
in some detail in later sections. As the properties of the anharmonic oscillator are well-known, we do not present 



simulations of the one-mode equations (14), but proceed directly to the multi-mode field problem. 



D. Multi-mode problem 

By analogy with the single-mode problem, in which single-mode operators were replaced by classical variables 
driven by white noise sources, we might expect that a complete quantum field can be replaced by classical fields 
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suffering independent noise sources at every point in space. There are a number of ways of proving this claim. A 
straight-forward (if notationally cumbersome,) method is to expand the field in a complete set of modes and mode 
operators and proceed by a direct generalization of the method in the previous section [4^. However, a more concise 
derivation is obtained by introducing the functional P-distribution |39|] 

P({V^,V*},t) = p('^)({V',^n,^)| (15) 

— 

where denotes the density operator p(t) antinormally ordered with respect to the field-operators in the 

Schrodinger picture. Putting the master equation obtained from the Hamiltonian (^) into antinormal order, and using 
the following functional analogues of the operator correspondences ^ : 

P-^ ^ ( V' - 77T ) -PW, p-^t ^ ■0+P(V'). (16) 



one finds the functional Fokker-Planck equation, 



As anticipated by the results for the single mode problem in section p C| , the diffusion matrix of this equation is non 
positive-definite and so there is no straight forward mapping onto a single stochastic differential equation ||l^. Just 
as for the single mode case, we move to a positive-P representation and double the phase space with the mapping 

ip{x,t) ipi{x,t), 
^+{x,t)^Mx,t), (18) 

where ipi{x,t) and rp2ix,t) are independent fields. As before, we are guaranteed a positive-definite diffusion matrix 
and finally obtain the pair of Ito stochastic equations 



idri'iix, t) — JC^i{x, t) + T^2{x, t)')PI{x, t) + ^/i^'^l{x, r)?7i(x, t), (19a) 
idrip2{x,T) = -)C'ip2{x,T) - T^pl{x , T)'ipl{x , t) + \/-ir?/'2 (a;, t)?72 (a;,T), (19b) 



where the noise sources 771 and 772 are real, Gaussian and delta-correlated in time and space: rji^x, T)r]j {x', r') = Sij5{x— 

x')8{t — t'). Note that the mean number of atoms {N{t)) = {f^^dxt[;^{x,T)ifi{x,T)) = ^'2(2;, t)V'i (a;, r) is 

conserved in the ensemble average but fluctuates during a single trajectory due to the complex noise. We remark that 
in practice, it is numerically more convenient to work with the complex conjugate equation to Eq. ( |19h| ). We have 
used the form shown in order to make a clearer connection to Eqs. (|l^). 

Although our derivation indicates that Eqs (^9|) allow one to calculate single-time normally- ordered quantum field 
averages, it was shown by Drummond |45| that they actually allow for multi-time time-normally- ordered averages to 
be found. In general, an expression for an arbitrary time-normally-ordered average is obtained by replacing 'ijj{x,t) 
by ipi{x,t), Tp^^x^t) by 'ip2X,t), and the quantum averaging by the stochastic one: 

(fV;t(a;,i)-->^(a;',t') TV>(a;", t") ■ ■ ■ ^{x'" , t'")) = V2(x, t) • • • ^-2(2:', t')Mx" , t") ■ ■ ■ M^'", *"')• (20) 

Here, T and T denote, respectively, direct and reverse time ordering of the field operators. The upper bar on the RHS 
of this relation denotes an averaging over the random trajectories {ipi, 1P2}, with the stochastic measure characterised 
constructively by Eqs (p^). In other words, this is a path integral over trajectories {tpi,ip2}, with Eqs (|l9|) providing 
the measure over the paths. In quantum field theory, quantum averages of the form in Eq. ( pO[ ) appear in the well- 
known Keldysh diagram techniques j46| . They are a subset of the full set of Keldysh averages (which in general also 
contain ip^s under the T-ordering and t/i^'s under the T-ordering). For this subset of quantum averages, Eqs ( p^ are 
fully equivalent not just to the master equation, but to the Heisenberg equations of motion for the quantum field, 
providing a constructive path-integral representations for these averages. (Moreover with external sources added to 
Eqs (|9|) they account for the full set of Keldysh averages, thus becoming fully equivalent to the Heisenberg equations. 
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see [^8[^). In the context of field theory, it is perhaps worth remarking on a helpful simplification that results from 
the non-relativistic nature of the problem. The density matrix at a time t — tq may be mapped directly onto the 
P, Wigner or positive-P distributions at the same time. These are then used as distributions for initial conditions in 
simulations. In the general case, one must rather match the density matrix to the distributions at r = —oo, subject 
to the usual assumption of adiabatic turning on of the interaction | |3^ . 

The positive-P representation is guaranteed to give exact results for as long as the ensemble averages converge. 
However, we see below in section |^ that the trajectories are prone to large excursions from the mean which quickly 
cause the simulation to blow up. Such problems with the positive-P representation are well-known and 
occur especially in systems with strong nonlinearity and weak (or vanishing) damping which is precisely the situation 
in the present case of a trapped interacting Bose condensate. We believe this is the first case however, for which 
divergent trajectories appear for realistic physical values. As such it is indication of the likelihood of strongly non- 
classical behavior outside the description of the GPE. It is important to realize that the failure of the positive-P in 
such cases is not indicative of a genuine "divergence" in the sense of quantum electrodynamics but merely represents 
a rapid (presumably exponential) growth in the width of the distribution. So while in theory the distribution remains 
physically correct, in practice it becomes impossible to accurately sample the whole distribution numerically. In fact, 
this problem is strongly dependent on the parameter range chosen. Drummond and Corney have successfully used 
the positive-P representation to simulate evaporative cooling of the condensate . The cooling problem is a case 
in which the positive-P representation can be expected to be more robust than in zero-temperature calculations for 
two reasons: the atomic density is much lower, at least initially, and there is a considerable damping in the form of 
the RF field used to remove the hotter atoms. 



E. Truncated Wigner representation 

In the absence of exact stable methods, we are forced to consider approximate simulation techniques. One approach 
that h as pro ved successful in optical problems is the "truncated Wigner" method [Q. In a similar fashion to that of 
section II D, the master equation can be mapped onto a Fokker-Planck like equation for the Wigner distribution | p^ , p^ , 
which returns symmetrized expectation values as opposed to the normally ordered averages of the P representations. 
That is to say, we define the Wigner distribution by analogy with Eq. (^5|) as 

Wiiij, r}, r) = p(^''"^-)({Vi, ^^}, r)| (21) 

where p^^^™-^ denotes the density operator p{t) symmetrically ordered with respect to the field-operators . (see 
Ref. [ [14[ for the connection of this definition to more familiar expressions for the Wigner function) . Using the functional 
differentiation notation, the operator correspondences take the form |Q, 

p^^{^-\^)w, f'^' - (22) 

Using these relations in the master equation we find the Wigner function evolution equation 



In this case, there is no second derivative term — the diffusion matrix vanishes identically — and the quantum noise 
acts via third order derivatives as "cubic noise" . Unfortunately, there is no simple mapping from cubic noise to a 
stochastic representation |l^ and as we have discussed earlier, a direct integration of Eq. ( p3| ) is impractical. The 
simplest approximation is to truncate Eq. ( p3| ) at second order so that we are left with a single deterministic equation 
for the classical field V'VKi which is just the standard time-dependent GPE 

idripwix^r) =K.iPw{x,t) +T\'ipwix,T)\'^il^wix,T). (24) 

Although this equation is completely deterministic, it is not the case that we have discarded all effects of quantum 
noise. Noise is still included explicitly in the initial state which is now represented as a distribution of functions 
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tpw{x, 0). (In fact, even in the positive-P representation we would require a random distribution of starting functions 
unless the initial state was a coherent state.) We discuss the choice and representation of initial states in detail in 
section pl| 

In optical problems, it has been found that the Wigncr approach gives accurate results in the large photon number 



limit, when it might be expected that the influence of the third order quantum noise is small. In section IV, we 
test the Wigner predictions against the positive-P results for as long as the latter are stable. Using the Wigner 
distribution entails one further limitation. Typically, the physically most interesting quantities are time-ordered, 
normally-ordered averages, as provided directly by the P representations ||l^ , [l5|| . As the Wigner distribution returns 
symmetrized moments and we do not know the unequal time commutators for the field operators, we can not usually 
find multi-time averages with the Wigner method. An exception is the two-time normally-ordered correlation function 
for coherent initial states, 

{^^x, r)^ix\ t')) = (^op(a;, r)| i^\x, t)^{x' , r') |^op(x, t)) 

= i/'Gp(a;,r)* (?/'Gp(a;,T)| ?/'(a:',r') \%I)c_i-{x,t)) 



= i}}q^{x,t)*'iPw{x,t'), (25) 

which is thus reduced to a single time expectation value with no ordering problems. Note that even for coherent 
initial states, higher-order correlations such as (: ?/'^(ti)'(/'(ti)V'^('^2)V'('''2) :) are unavailable. 



III. INITIAL STATES 



The question of suitable initial states for simulation is somewhat involved. Here we wish to use states that can 
be thought of as a good representation of the "ground state of the condensate", in as much as this is possible in 
a symmetry broken picture. We consider only zero temperature states here. For T = 0, the simplest option is to 
choose an initial coherent state, that is, precisely that state assumed at all times in a conventional calculation with 
the GPE. Our simulations then indicate how the actual state evolves away from the coherent state. To do so, we set 
the mean field {il]{x, 0)) equal to the solution of the time-independent GPE i!or{x) (the "ground state wavefunction" ) 
and assume vacuum noise in all modes. For the normally-ordered positive-P representation, vacuum noise is obtained 
simply by the choice tpi[x,Q) = 'ip2{x,0) — V'Gp(a^)- In the symmetrically-ordered Wigner representation, the noise 
must be explicitly included in each mode of a suitable basis. Each trajectory begins with a different field of the form 

N 1 

ijw{x,0) = tpGpix) + ^ ■^Vj(t>jix) (26) 
i=o 

where rjj is a complex random variable of zero mean with rjjrjk — and r]*rik — 5jk. The sum is taken over N modes 
of a complete basis {^^(a;)} and N is taken sufficiently large that the results are independent of N. Natural choices 
for the {^fc(x)} are the discrete position basis (pjix) = 5{x — jAx)/y/'Ax or the harmonic oscillator basis. The latter 
has the advantage that if N is not too large, the modes do not extend to the boundary of the simulation window 
and there is no risk of noise artificially "wrapping around" the simulation. In practice, we have seen no difference in 
results when using either of these bases. 

However, if one wishes to find a good approximation to the ground state of the many body system, the coherent state 
is certainly not an optimal choice. Of course, as we are using a symmetry-breaking approach, no state can be truly 
stationary — there must always be a degree of phase diffusion associated with the number superposition implied by the 
assumption of a non-zero mean field. Nevertheless, there are states we might favour over the coherent state. Standard 
applications of Bogoliubov theory at zero temperature [[l9[|5C| ] approximate the second quantized Hamiltonian by the 
diagonal expression 

H = K + Y,Eib]bj, (27) 

where K is a, constant, hj is the annihilation operator for the quasiparticle excitation of energy Ej and the mean field 
satisfies the GPE. Further, the field operator may be written as 

i^{x)=Y.^u,{x)b,~v*{x)h% (28) 
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where the mode functions Uj and Vj are solutions to the BogoUubov-de Gennes eigenvalue equations 1^,^ . Our 
second choice for the ground state is thus the vacuum in the Bogoliubov representation: 

N 1 

il^w{x,(}) = V'gpW + ^ i^iVjUjix) - VjVjix)). (29) 

However, Lewenstein and You [ pT| have pointed out that in a symmetry-breaking Bogoliubov method, the existence 
of a zero-energy Goldstone mode requires the inclusion of an extra term in the Hamiltonian involving the condensate 
"momentum" P that accounts for the phase diffusion of the mean field. In this case we have 

H = K+^P'+J2Eibyb,, (30) 

j 

with a = Ndfi/dN, P = J^^dx4>Gp{x){5'ip + Sip^), and S'lp = tp — {^). This Hamiltonian implies an infinite 
amplitude squeezing of the condensate which is clearly unphysical. It has been shown elsewhere for one and two 
mode models |5^,^ that retaining cubic and quartic terms in the Hamiltonian which are neglected in the Bogoliubov 
method leads to a finite squeezing. Here, we do not perform a full treatment of the effect of the higher order terms 
for the multi-mode system, but take for our third "ground state" , the lowest-energy variational state in which each 
mode in the Bogoliubov basis is independently in a minimum uncertainty Gaussian state. 

We also briefly remark that the choice of initial state is closely tied to the manner of state preparation. In many 
instances, the appropriate state need not be the ground state. In a recent experiment at JILA a single condensate 
is subjected to a short tt/2 pulse creating a second condensate in a different internal state. As the pulse length is 
shorter than the time required for significant nonlinear dynamical effects to occur, the combined two condensate 
system might be expected to exhibit binomial statistics. If the trapping potential were arranged so that the two 
clouds did not subsequently overlap, we could then model the evolution of one of the two condensates assuming a 
number variance (AiV)^ « N/2. In fact, it may be checked using results in Ref. |^l[, that the state in which all the 
Bogoliubov modes are in the vacuum has number statistics for the condensate mode very close to (AA^)^ « N/2. 
For a slower transfer of population, nonlinear effects would play a role and a more complex initial state would be 
appropriate [Q. 

As an example of a completely different initial condition, the first-principles simulation of evaporative cooling using 
the positive-P representation that was mentioned earlier |4^, begins essentially with a thermal state for the atom 
field. 



IV. RESULTS 
A. Numerical methods 

The parameters for our simulations are chosen to represent the following system. We consider a condensate of 
N — 1000 sodium atoms in a cylindrical trap with A = w^/w^ = 0-025 so that the one-dimensional approximation 
is reasonable. The radial frequency is set at either LOy./2'K = 800 Hz (the "strong trap") or uj^/2'k — 200 Hz ("weak 
trap"). Taking the scattering length as a = 4.9 nm we obtain for the nonlinear constant F = Fgtrong = 0.084 or 
F = Fwcak = 0.042 . The initial state mean field solutions were obtained by imaginary-time propagation of the GPE 
and quasiparticle energies and mode functions found by standard methods ^,|5^ . The predictions of our simulations 
were checked by ensuring that the results did not change when the time step was decreased or the size of the spatial 
grid increased. Simulations in the truncated Wigner representation were performed with a standard second-order 
split step method jij]. Due to the large nonlinearity of the system and consequently strong noise in the positive-P 
simulations, the standard Euler split-step algorithm was not able to give results independent of time step even at a 
step size of drr = 0.0005. Hence, we used a strongly-convergent semi-implicit method which gave reliable results 
with 256 spatial points and a time step dr = 0.001. 
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B. Quantities of Interest 



We present our results in terms of three general quantities. To demonstrate the ability to determine two-time 
correlation functions we would like to calculate the quantity 



9^'\r,0) = ^ L (31) 



This is straight-forward in the positive-P representation. However, as discussed in section [IE, in the Wigner rep- 
resentation we may only calculate unequal-time normally ordered correlations for coherent initial states, when from 
Eq. (p5|) we have 



gW(.,0)^-^-°°-^"^SpWM.(a.,r)^ (32) 

For other initial states, we can define a nominal "condensate mode" operator associated with the normalized solution 
to the GP equation, 

/oo 
dxij(.p{x)il>{x,t), (33) 
-oo 

where ijj^j,(x) = ipcrix)/^ J^^dx \'4'gp{x)\^ ■ Its mean value (aGp(i)) still monitors the collapse of the wave function 
but is strictly a one-time average and can be calculated in either representation.We may also calculate the occupation 

of the condensate mode (nop) = (aGpaop 



Finally, spatial correlations may be analyzed in terms of a spatial squeezing spectrum. We define the localized 
amplitude quadrature operator 

X0{x,t) = ^^^{x)i^^{x,T)e''' + ^^^{x)4,{x,T)e-'''. (34) 
Defining the Fourier transformed operator 

Xe{k,T)^^ dxe*'="Xe(x,T), (35) 



the squeezing spectrum is defined as the normally-ordered expression [p^ , p8|j2g| , |55|] 

Se{k,T) = 2t: Xe{-k)Xg{k) ■}j , (36) 
which in the Wigner representation becomes 



Se{k,T) = -l + 2T:Xe{-k)Xe{k). (37) 

The angle 9 for optimum squeezing is in general a function of k. Hence a useful quantity is the spectrum of "best 
squeezing" S-[aa_^{k^T) which gives the largest possible squeezing at each wavenumber component k. 



C. Comparison of methods 



We first give some examples of calculations with the same parameters using both the positive-P and Wigner 
simulations. Figure ^ shows the two time correlation function g^^^ (r, 0) for 1000 atoms in the strong trap configuration 
for a coherent initial state. Single mode models predict a Gaussian decay 

g^^\T, 0) = cxp[-aV2(A7V)V(2Ar2)] (38) 

where a — Nd^/dN, and (AA^)^ = iV is the variance in atom number of the initial state. This model is indicated 
by the chain curve. The Wigner prediction (dotted) roughly follows the Gaussian decay but shows slow oscillations 
about the single mode curve, and in particular exhibits a linear decay at short times. The Wigner method gives a 
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stable result for arbitrary times. In the inset we show the short time behaviour, with the inclusion of the positive-P 
prediction in the solid line. This line stops at just r = 0.3 at which point unstable trajectories appeared. Also just 
visible are error bars on the positive-P line denoting one standard deviation (OSD) uncertainties. The dashed curves 
indicate the OSD errors for the Wigner calculation. The two methods clearly agree up to the point at which the 
unstable trajectories arise. Note that the positive-P error bars are very small right up till that point, indicating the 
sudden rapidity with which the distribution diverges. In Fig. ^ we show the occupation of the GP mode ( 
function of time with the Wigner result shown as the dotted line and the positive-P result shown as the solid line. 
We see oscillations in the number with an amplitude of around 5 % of the starting population. Once more while the 
Wigner result is stable, the positive-P simulations fail in a very short time. Note also that, the blow-up occurs very 
rapidly — there is very good agreement until just before the fatal moment with OSD errors for both methods being 
smaller than the thickness of the lines. 

Thus for the trap parameters considered so far, the positive-P representation is effectively useless. While instabilities 
of the positive-P are well-known this is perhaps the first occasion in which they arise in experimentally accessible 
parameter ranges. Given that the positive-P fails well before the completion of a single oscillation in Fig. ||, we 
might wonder how closely the truncated Wigner results approximate the true dynamics. One approach is to perform 
simulations at artificially low scattering lengths for which the nonlinearity is less severe and the positive-P simulations 
more robust. In fact the observation of Feshbach resonances in an optically trapped Na^'^ condensate fs^ demonstrates 
that reduced scattering lengths are now attainable in the presence of a sufficiently strong magnetic field. Figure ^ 
shows the GP mode occupation as a function of time for 1000 atoms in a coherent initial state, with the reduced 
interaction F = Pstrong/lO. The line styles have the same meaning as in Fig. |2| The positive-P trajectories are now 
stable for much longer and it is seen that both methods produce oscillations that are in agreement within the error 
limits. Note that the error limits grow in time for the positive-P but remain approximately constant for the Wigner, 
corresponding to the fact that no new noise is added after the initial condition for the Wigner method. We can thus 
now have some confidence that the Wigner calculations give results that are reasonably accurate for relatively large 
condensates. 



D. Comparison of initial states 



We now examine the behavior exhibited by different initial states. As explained in section [II, we compare the 
standard choice of coherent initial state with the vacuum state in the Bogoliubov representation and the Gaussian 
state with independent squeezing in each Bogoliubov mode. Figure Q shows the mean amplitude in the GP mode 
(ogp (■'")) with Wigner results shown in the solid lines and the single mode estimates based on the initial number 
variance shown in the chain lines. The mean amplitude is apparently described relatively well by the single mode 
model. The differences between the curves is largely accounted for by the difference in number variance in the three 
cases, which had the values {AN)'^/N — 1, 0.5 and 0.12 in the coherent state, Bogoliubov vacuum and squeezed 
Bogoliubov vacuum cases respectively. Figures give the spectra of best squeezing for the three initial states 
plotted as the function log[l -I- S'niax(fc, t)]- For the coherent state in Fig. ||, there is of course initially no squeezing. 
For a short time, there is significant squeezing at low wavenumbers. However, at large times the phase diffusion causes 
the long wavelength fluctuations to grow without limit jSlj l and the squeezing is destroyed. The other two initial states 
shown in Figs, ^and ^ show similar trends at large r but are clearly different at early times when the statistics of 
the initial state have not been swamped by the phase diffusion. This suggests that the squeezing spectrum may be 
a useful way of characterizing different quantum states of the condensate. Note that a single mode model could not 
predict different rates of change for the squeezing at different wavenumbers. 



E. Negative scattering lengths 



Finally we briefly examine the dynamics for a single case with a negative scattering length. In this case, the 
attraction between the atoms leads to a high density at the center of the trap, and consequently the nonlinear terms 
play a stronger role than in the positive scattering length case. To avoid the need of extremely fine spatial and 
temporal grids, we therefore use the parameters of the weak trap, with the scattering length set at Oneg = "QNa/lO 
giving r = —0.0042. The two-time correlation function ^'^'(t, 0) and occupation (ncp) display similar behavior to 
that seen earlier for the positive scattering length. Here we concentrate on the squeezing spectrum which is shown in 
Fig. Ij. This figure shows quite different structure to the earlier squeezing spectra with strong anti-squeezing at for 
wavenumbers near fc « 3, corresponding to a length scale of the condensate or "soliton" width. In fact, this spectrum 
is very similar to the spectrum of best squeezing for a fiber soliton with a Kerr law nonlinearity pqj29|| . This is not 
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surprising. With a strong negative nonlinearity in a one-dimensional trap, the condensate becomes strongly localized 
at the bottom of the trap. The nonlinearity dominates over the trapping potential and the ground state wavefunction 
is well-approximated by the fiber soliton expression V'(a^) = sech(V NTx), with a slight additional confinement due 
to the potential. Then as the propagation equations for the two systems differ only by the inclusion of the potential 
for the condensate, we can expect virtually identical spectra. 



V. CONCLUSION 



In this paper, we have applied phase space techniques for the propagation of a complete quantum field to the 
problem of a one-dimensional trapped Bose-Einstein condensate. As such systems are highly nonlinear and weakly 
damped, the exact approach using the positive-P representation is useful only for short times compared to the trap 
period and we are forced to use the approximate truncated Wigner method. For parameter ranges in which both 
methods work, we find agreement between the two. The Wigner method is stable and allows the calculation of 
one-time averages and certain conditional multi-time averages over long periods. Dynamics may be calculated for 
virtually any initial state with a reasonably well-localized Wigner function. 

It is interesting to compare our approach here with another set of tools for discussing quantum statistical properties 
of condensates — the rapidly growing field of Quantum Kinetic Theory (QKT), which has been developed in particular 
by Gardiner and ZoUer and their co-workers [^,^. In QKT, the system is divided into two distinct parts — the 
"condensate region" consisting of the condensate itself and a considerable number of the low-lying excitations and 
the "thermal region" which is essentially everything else and acts as a reservoir for the condensate region. One 
may then obtain master equations for the condensate region of varying complexity based on assumptions about the 
exchange of atoms between the condensate and reservoir. In our own approach, there is no distinction at all into 
condensate and thermal atoms and thus no approximations required in order to implement such a distinction. The 
condensate itself plays no privileged role within the model and we work simply with one complete quantum field. The 
special properties normally associated with condensates are manifested just as different correlations of the quantum 
field. The stochastic method described here thus may also serve to provide comparisons with the predictions of QKT 
from a rather different vantage point. Indeed, Drummond and Corney's simulations of evaporative cooling using the 
positive-P representation have produced |Q similar results to kinetic models of evaporative cooling |58|. Moreover, 
one might envisage a hybrid model in which the low-lying condensate modes are treated using a stochastic approach 
while the upper modes are reduced to a thermal reservoir using the techniques of QKT. 

Finally, we point out some of the systems to which this theory could be easily applied. As mentioned earlier, the 
coherence properties of the output beams of atom lasers are certain to be of central importance in the near future. 
Calculation of two-time correlations and squeezing spectra for various laser designs is a natural application. The 
phase diffusion between coupled condensates is also beginning to attract interest and has currently only being studied 
theoretically within the context of Bogoliubov theory 26 59[. 
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FIG. 1. Two-time correlation function g^-^^r) for coherent initial state in the strong trap. Line styles are positive-P (solid), 
Wigner (dotted), single-mode model (chain). Inset shows early times with 1 standard deviation errors for the Wigner method 
shown in dashed. 
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FIG. 2. Occupation number (ncp) as a function of time for strong trap. Line styles are positive-P (solid) and Wigner 
(dotted). 



FIG. 3. Occupation number {hop)- Line styles are positive-P (solid) and Wigner (dotted). 



FIG. 4. Mean value (aGp(t)). Wigner results are shown in solid lines, 1 mode models in chain. 



FIG. 5. Maximum squeezing spectrum plotted as log[l + iSmax] for weak trap parameters with coherent initial state. 



FIG. 6. Maximum squeezing spectrum plotted as log[l + Smax] for weak trap parameters with Bogoliubov vacuum initial 
state. 



FIG. 7. Maximum squeezing spectrum plotted as log[l-|-ySmax] for weak trap parameters with squeezed Bogoliubov vacuum 
initial state. 



FIG. 8. Maximum squeezing spectrum S'max for system with attractive interactions: scattering length Oneg = —0.049 nm 
and N = 1000 atoms. 
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